Broken-symmetry v — quantum Hall states in bilayer graphene: 
Landau level mixing and dynamical screening 
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For bilayer graphene in a magnetic field at the neutral point, we derive and solve a full set of gap 
equations including all Landau levels and taking into account the dynamically screened Coulomb 
interaction. Lhere are two types of the solutions for the filling factor v = 0: (i) a spin-polarized 
type solution, which is the ground state at small values of perpendicular electric field Ej_, and (ii) 
a layer-polarized solution, which is the ground state at large values of E± . The critical value of E± 
that determines the transition point is a linear function of the magnetic field, i.e., E± jCr — _E° ff + aB, 
where _E° ff is the offset electric field and a is the slope. The offset electric field and energy gaps 
substantially increase with the inclusion of dynamical screening compared to the case of static 
screening. The obtained values for the offset and the energy gaps are comparable with experimental 
ones. The interaction with dynamical screening can be strong enough for reordering the levels in 
the quasiparticle spectrum (the n — 2 Landau level sinks below the n = and n = 1 ones). 

PACS numbers: 73.22.Pr. 71.70.Di, 71.70.-d 
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I. INTRODUCTION 

Bilayer graphene is a unique material in con- 
densed matter physics. It combines some characteris- 
tics of monolayer graphene and more traditional two- 
dimensional electron systems. Its low energy electron 
spectrum is gapless and is given by parabolic valence 
and conduction bands with massive chiral charge carri- 
ers touching at two K and K' valley points. As was first 
suggested in Ref. [l] and experimentally shown in Ref. 2, 
if a potential difference between the layers is applied, a 
tunable energy gap is opened at the valley points. In 
fact, as was indicated in Ref. 0, even without an external 
electric field, the quadratic dispersion relation in bilayer 
graphene implies that the electron-electron interaction 
should open a gap in the spectrum at the neutral point 
in clean bilayer samples (for a similar development in the 
theory of topological insulators, see Ref. 0). 

The rich spin- valley approximate SU (4) symmetry of 
the low energy electron Hamiltonian with the Coulomb 
interaction leads to many interesting possibilities for the 
gap generation. In the absence of magnetic field, anoma- 
lous quantum Hall (QAH)^£ quantum spin Hall (QSH), 
layer antiferromagnet (LAF), and nematio 7 -^ states were 
suggested as possible ground states of bilayer graphene 
at the neutral point (for a general discussion, see Ref. @) . 

As is well known, a magnetic field is a strong catalyst 
of symmetry breaking in graphene like systems J^r— In 
the presence of a magnetic field, the gap generation was 
experimentally observed in Refs. ll3l - f20L It was found 
that the eightfold degeneracy in the zero-energy Landau 
level can be lifted completely, giving rise to the quan- 
tum Hall states with filling factors v — 0, ±1, ±2, ±3. In 
suspended bilayer graphene used in Refs. llWlaflol the 



values of the gaps are of the order of a few milli elec- 
tron volt for magnetic fields B ~ 1 T. The theory of the 
quantum Hall (QH) effect in bilayer graphene has been 
studied in Refs. l2ll - [3ll . 

It was revealed in Ref. [TH that the energy gaps scale 
linearly with magnetic field B in bilayer graphene. This 
is in contrast to the case of monolayer graphene where 
a \f~B scaling for the gaps takes placed As was sug- 
gested in Refs. I25U27I a strong screening of the Coulomb 
interaction is responsible for this modification of the 
scaling in bilayer. The physics underlying this effect 
is the following ! 26 : 27 The polarization function in bi- 
layer graphene is enhanced much more than in monolayer 
graphene and, as a result, its contribution to the effec- 
tive interaction dominates over that of the bare Coulomb 
potential for mag netic fields B < 30 Such a strong 
screening radically changes the form of the interaction 
and leads to a linear scaling law. 

Another interesting phenomenon in the v = QH state 
in bilayer graphene is the phase transition between the 
spin polarized (ferromagnetic) phase and the layer po- 
larized one in the B-E± plane, where E± is an electric 
field orthogonal to the bilayerplanes. It was predicted 
in theoretical studies in Refs. [26H29I and observed in ex- 
periments in Refs. [l5l fl7lfl9l 

In suspended bilayer graphene used in Refs. [l3l . 
[lsIIlgI . the mobility is in the range from 10, 000 to 
15, 000cm 2 /Vs and the values of the gaps are of the or- 
der of lmeV for magnetic fields B ~ IT. Note that 
the energy separation between the lowest Landau levels 
(11 = 0, 1) and the n = 2 level in the free bilayer model 
in a magnetic field is ^/2Tiuj c ^ where the cyclotron en- 



ergy is hui c = h\eB\/mc c 
m of quasiparticlcs is m 



2.155 [T] meV and the mass 
~ 0.054m e . Even though it 
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is of the same order of magnitude as the experimental 
gaps, the energy separation between the Landau levels 
may be large enough to justify the use of the lowest Lan- 
dau level (LLL) approximation in the analysis of the gap 
equations i 2 1 ' 25 ~ J2£ 

Recently, however, the experiment a 18 ^ 19 in suspended 
bilayer graphene with a higher mobility (80, 000 to 
100, 000 cm 2 /Vs) revealed much larger gaps, about 
6.5 me V at B ~ IT. It is natural that a higher mo- 
bility specimen produces much larger gaps compared to 
those in Refs. Il3lll5lll6l . Formally, such large gaps exceed 
the energy of the 3rd Landau level in the free low-energy 
effective theory of bilayer grapheneji i.e., in a model with 
the Coulomb interaction ignored. Obviously, in this case 
the use of the LLL approximation in the study of the 
generation of the gaps cannot be justified. In this study, 
therefore, we amend the analysis of Refs. I26ll27l by prop- 
erly taking Landau level mixing effects into account.— 

The inclusion of Landau level mixing is also instruc- 
tive from the viewpoint of comparing the quantum Hall 
ferromagnetism2i (QHF) and magnetic catalysis^ (MC) 
scenarios describing the gap generation in graphene (for 
example, see discussions in Refs. I36I439I ) . Technically, 
the difference between them comes from the choice of or- 
der parameters describing the SU (4) symmetry breaking. 
While the QHF and MC order parameters are inequiv- 
alent in general, they are indistinguishable in the LLL 
approximation ! 2 ^ 37 ! 38 When the gaps are large and Lan- 
dau level mixing is essential, one cannot avoid/neglect 
any of these order parameters. Therefore, by solving the 
corresponding gap equations and determining the ground 
state, one has an opportunity to shed some light on the 
quantitative roles of the QHF and MC in the dynamics 
of graphene. 

It is interesting to compare the importance of Landau 
level mixing in bilayer graphene with that in monolayer 
one.— ~— . The study in Ref. revealed several qualita- 
tive effects in the latter, including the "running" of the 
gaps and the Fermi velocity with the Landau level in- 
dex n. A detailed analysis of the quasiparticle spectra 
in higher Landau levels (n > 1) also reveals the role of 
QHF and MC order parameters. The same is expected 
in bilayer graphene. However, the Coulomb interaction 
in bilayer should play a more profound role. This is the 
consequence of a much smaller characteristic energy scale 
hLij c in bilayer graphene as compared to the Landau en- 
ergy scale Et ~ 26 \JB [T] meV in monolayer one. 

In this study, we make another essential improvement 
in the theoretical analysi o 26 i 27 ' 31 by replacing the static 
(i.e., instantaneous) screening of the Coulomb interaction 
by the dynamical one. As will be shown below, these two 
new ingredients lead to both larger gaps and an offset 
for the critical value of the electric field £?° > which sep- 
arates the spin-polarized phase and the layer-polarized 
one at zero magnetic field. The latter is observed in all 



experiments in suspended bilayer graphene i 13 ' 15 i 16 ' 19 

Another interesting finding of our study of the dynam- 
ics with dynamically screened Coulomb interaction in the 
v = QH state is a Landau level reordering, caused by 
large gaps in the lowest Landau levels. The result is a 
dramatic rearrangement of the quasiparticle spectrum: 
the lowest conduction band and the highest valence one 
are now given by the n = 2 Landau level, rather than by 
the LLL (n = 0,1). 

This paper is organized as follows. In Sec. [H] 
we introduce an effective low-energy model of bilayer 
graphene. In this section, we also numerically calculate 
the frequency-dependent polarization function. We find 
that the corresponding function allows a simple fit given 
by the product of the static function and a simple form 
factor depending on the energy and momentum. The po- 
larization function is used in a coupled set of gap equa- 
tions derived in Sec. IIIII In Sec. IIV1 we present our nu- 
merical results. A general discussion of the main results 
and their comparison with experiment are presented in 
Sec. El Several appendixes at the end of the paper con- 
tain technical details and derivations used supplementing 
the presentation in the main text. 



II. MODEL 

We will utilize the same model as in Refs. I26ll3~i1 . The 
free part of the effective low-energy Hamiltonian of bi- 
layer graphene readsi 

where it = p x + ip y and the canonical momentum p = 
— iftV+eA/c includes the vector potential A correspond- 
ing to the external magnetic field B = (0||, B±), which is 
taken to be orthogonal to the bilayer planes, B = \B±\. 
The quasiparticle mass is m = ji/2vp f=a 0.054m e , where 
vf ~ 8.0 x 10 5 m/s is the Fermi velocity, 71 0.39 cV, 
and m e is the mass of the electron. The two compo- 
nent spinor field \&f s carries valley (£ = ± for valley K 
and K' , respectively) and spin (s =4-, t f° r spin down 
and up, respectively) indices. We will use the standard 
convention^ = (ipKA 1 ,'*pKB 2 )s f° r valley K and 

^-s = {i } K , B 2 i ipK'A^s for valley K' . Indices A\ and B2 
label the corresponding A and B sublattices in the layers 
1 (top) and 2 (bottom), respectively, which, according 
to Bernal (A2 — B\) stacking, are relevant for the low 
energy dynamics. Let us emphasize that the sublatticc 
and layer degrees of freedom arc not independent in this 
low-energy model: the sublattices A and B correspond 
to the layers 1 and 2, respectively. 



The Zeeman and Coulomb interactions plus a top-bottom gates voltage imbalance Aq in bilayer graphene are 
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described as follows: 



Za 3 + A £t 3 



Aq | 
77171 ' 



TT^IT 
— 7r7T 



H int =Y,J d 2 r^l(r) 

+ 1 1 d 2 rdV {v> - r') [ Pl (r) Pl (r') + p 2 (r)p 2 (r')] + 2^ 12 (r - r')pi(r)p 2 (r')} , 



(2) 



Here Z = p, B B = 0.027 hoj c = 0.059£[T] meV is the 
Zeeman energy, er 3 is the diagonal Pauli matrix in spin 
space, and r 3 is the diagonal Pauli matrix acting on the 
two components of the fields and ^- s - Note that the 
presence of £ in the voltage imbalance term is related to 
the different order of the Ai and i? 2 components in $! +s 
and "f-s [compare also with the layer projection opera- 
tors in Eqs. ((3j) and (|4]) below]. The voltage imbalance 
between the top and bottom gates is related to the elec- 
tric field applied perpendicularly to the bilayer planes: 
A = eE±d/2. 

The third term in the square brackets in the first line 
of Eq. © leads to a small splitting of the Landau levels 
with orbital indices 77 = and 77 = \X However, due 
to the factor 71 in the denominator, it is suppressed in 
comparison to other terms in the interaction Hamiltonian 
@. In fact, the splitting that it produces is of order of 
O.OlAo£>[T] that is substantially less than the splittings 
due to the other two terms in the square brackets in the 
first line of Hamiltonian © for reasonable values of a 
magnetic field B (for a recent discussion of this term, 
see Ref. l3l[ ). As will be shown below, it is also much 
less than the dynamical splitting between the 77 = and 



77 = 1 levels due to the Coulomb interaction. Because of 
that, it will be omitted in the analysis below. 

The Coulomb interaction term V(r) in H lnt is the 
bare intralayer potential whose Fourier transform is given 
by V(p) = 2-Ke 2 /(np), where K is the dielectric con- 
stant. The potential T4 2 (r) describes the interlayer elec- 
tron interactions. Its Fourier transform is Vviijp) = 
2ire 2 e~ pd / (np) : where d = 3.5 x 10~ 10 m is the distance 
between the layers. The two-dimensional charge densities 
in the two layers are 

Pl (r) = X>S.( r )*MO**.(r), (3) 
P2 (r) = X>S.( r )*MO<Mr), (4) 

where 7>i(0 = C 1 + ^ 3 )/ 2 and ^MO = (1 - £t 3 )/2 arc 
projectors on the states in the layers 1 and 2, respec- 
tively. When the dynamical screening effects are taken 
into account, the potentials V(r) and Vi 2 (r) are replaced 
by effective interactions V e g(t,r) and Vi2eff(i, r ) which 
are no longer instantaneous. 



The Schwinger-Dyson (gap) equation for the quasiparticle Green's function (propagator) in the Hartree-Fock ap- 
proximation reads^i 

G- 1 (t-t';r,r') = S -1 (i - t';r,r') - iG(t - t';r,r')V eS (t - t';r - r') 

- * [Pi(®G(t ~ t'; r, r')V 2 (0 + P 2 (flG(t - t'; r, r')Pi(f)] V lh (t -t';r- r') 

- \[Vi{0~V 2 {0] tT{[V 1 (O-V 2 (O]G(0;0)}V^(0)6(t-t')S(r-r'), (5) 



where Vi h (t -t';r- r') = V 12cS (t -t';r- r') - V cS (t - 
t'\ r— r'), S(t—t'; r, r') and G(t—t'; r, r') are the free and 
full Green's functions, respectively. These two Green's 
functions are described in Appendix [A] In the presence 
of an external magnetic field, they are not translationally 
invariant functions but can be written in the form of 
a product of a universal Schwingcr phase (which spoils 
the translational invariance) and translationally invariant 
functions. 



The momentum space expressions for the interaction 
potentials in the momentum space ar o 27 ' 43 



V?L arc (0) 



27re 2 



1 



« p+^fn(u :P y 



2ire 2 d 



(6) 
(7) 



The explicit form of the interlayer potential Vil(oj,p) can 
be found in Appendix A of Ref. [27]. However, here we do 
not need it: due to the presence of the projectors Vi(£) 
and 7*2{C) m the second line in Eq. ([5]), the corresponding 
Fock term does not contribute to the final form of the gap 
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equation. Note also that since bilayer has no net electric 
charge at the neutrality point, we dropped all the Hartree 
terms proportional to tr[G(0; 0)] (i.e., the total density of 
charge carriers) in the gap equation^ 

The screened Coulomb interaction V e g(uj,p) depends 
on the dynamical polarization function H(uj,p). In the 
previous theoretical studies of the gap equation in bilayer 
graphene in a magnetic field, only a static polarization 
function was used<2£ r 28 ' 31 In the random phase approx- 
imation, the static polarization function in a magnetic 
field was calculated in Rcfs. [2611271 and is rcplotted in the 
left panel in Fig. [T] The corresponding approximation 
for the interaction potential reads 



V cS (0;y) 



2ir£ 2 huj r 



Ky/xy + 4tiT%) 



(8) 



where y = (p£) 2 /2 is a dimensionlcss variable used in- 
stead of the wave vector p and 11(0, p) = (m/h 2 )Ii(y). 
By definition, £ = ^Jhc/\eB\ is the magnetic length and 
x = 2h 4 /(e 4 m 2 £ 2 ) w 0.0035 [T] is a dimensionless pa- 
rameter which determines the value of y below which the 
screening effects are negligible. 



The use of static screening approximation greatly sim- 
plifies calculations. Our analysis below shows, however, 
that this approximation significantly underestimates the 
strength of the Coulomb interaction in the Fock term. In 
fact, we find that taking into account the effects of dy- 
namical screening effectively leads to about three times 
larger gaps compared to the case of the static screening. 



In the one loop approximation, the frequency-dependent polarization function in Euclidean space (after the Wick 
rotation) is given byS 7 - 



nfoy) 



Mr 



irt^Ml + a 2 



[Co,n{v) +&,n(v)] 



M n + M r 



2lT ^ ( Mn + Mm) 2 +cT 2 

n.m—2 v ' 



(y) + £ n -2,m-2(y) 



M n M m 



,(2) 

-n-2,m-2 



(y) 



(9) 



where the dimensionless parameters M n = \J n(n — 1), a = oj/uj Ci and the following functions were introduced: 

41 = 2^/ d2re ~ tpr (J) (^) L l (^2) = j dtt a e- t .M2Wt)L r l(y)L: i (y) 



( _ ir+m (!!L^ e -, L „-™ (y)i m- n 



(10) 



Here L" (y) are the generalized Laguerre polynomials and 

Jo(x) is the Bessel function (£ n , m = Cn}m)- (For calcu- 
lation of such integrals, see Appendix A in Ref. H3-) In 
the case of a = 0, 2, in particular, one obtains 



£ n ,m(y) = (-1) 



n+m —y r n—m 



(y)L™- n (y), (11) 



As expected, the dynamical effects reduce screening at 
nonzero ui. 

We find that, in a wide range of momenta and frequen- 
cies, < {pi) 2 /2 < 10 and < u/uj c < 10, the dynam- 
ical polarization function can be well approximated (to 
within a few percent) by the following fit: 



4 2 L(y) = (-l) n+m e-y(m+l)(m+2)L^S (v) K~ n GO ■ 

(12) 

The above expression for the dynamical polarization 
function can be calculated numerically with a relative 
ease by employing a simple trick. First, we subtract the 
static part from Eq. ([9]) to obtain an expression for the 
difference n(er, y) — 11(0, y). Such a difference is given in 
terms of a rather quickly convergent series and, there- 
fore, can be tabulated as a function of two variables, a 
and y, without much effort. After adding the static part, 
which is a function of only one variable and which was 
calculated earlier ; 26 ' 27 we finally obtain the polarization 
function II (er, y). It is plotted in the right panel in Fig.[T] 



n fit (<r,2/) = n(o, y ) 



l + b lV 



y/1 + b 2 y + b 3 y 2 + b 4 a 2 ' 



(13) 



where the fit parameters are bi — 0.608, &2 = 1.572, 
63 = 0.357, and 64 = 0.868. Below this fit is used in 
our analysis of the dynamics to model the interaction 
potential with dynamical screening, i.e., 



2ir£ 2 huj r 



Ky/xy + 47rn fit (cr, y) 



(14) 



which will replace the approximation with static screen- 
ing in Eq. ©. 



5 




FIG. 1: (Color online) The static polarization function (left) and frequency-dependent polarization function (right). 



III. GAP EQUATION 

As discussed in the Introduction, here we consider a 
rather general ansatz for the full Green's function that 
includes the effects of both QHF and MC order parame- 
ters. While the former describe chemical potentials con- 
nected with conserved charges, the latter describe Dirac 
like masses. 

As discussed at length in Ref. H3, if both the Zee- 
man and Ao terms are ignored, the Hamiltonian H = 
Ho + H int , with H and Hi nt in Eqs. JT]) and @, possesses 
the symmetry G = U^ K \2) S x U^ K '\2) S x Z$ x Z$, 
where t/ (y) ( 2 )s defines the U(2) spin transformations in 
a fixed valley V = K,K', and Z^J describes the val- 
ley transformation £ — > — £ for a fixed spin (s =^,^). 
The Zeeman interaction lowers this symmetry down to 
G 2 = C/ (K) (l)t x U^ K \l) i x U^ K 'Hl) t x U^ K '\l) l x 

Z 2V x Z 2V' wnc rc U {v) (\) s is the U(l) transformation 
for fixed values of both valley and spin. Including the 
Ao term lowers the Gi symmetry further down to the 
G 2 = tf<*>(l)t x uW^x u( K )(l) t x U^ K '\l) { . 

The dynamics in the integer QH effect in bilayer 
graphene is intimately connected with dynamical break- 
down of the G and G2 symmetries. Two sets of the order 
parameters describing their breakdown were considered 
in Refs. [26ll27l . The first set consists of the quantum Hall 
ferromagnetism (QHF) order parameters^ 

& : E<*L^> ■ ( 16 ) 

s 

While the order parameter (fTo)) is a conventional ferro- 
magnetic one, the order parameter (|16p determines the 
charge-density imbalance between the two valleys. The 



corresponding chemical potentials are ^3 and /2 S , respec- 
tively. 

The second set consists of the magnetic catalysis (MC) 
order parameters^ i.e., the Dirac A s and Haldane A s 
mass terms: 

A s : £<4 s r 3 * ? , s ), (17) 

A < : E<*U r3 *^> ■ ( 18 ) 

The order parameter (|17[) describes a charge density wave 
in both the K and K' valleys. This order parameter pre- 
serves the G2 symmetry^ On the other hand, the order 
parameter (|18p . connected with the conventional Dirac 
mass A s , determines the charge-density imbalance be- 
tween the two layers^ The structure of this mass term 
coincides with that of the bare voltage imbalance term 
between the top and bottom gates introduced in Hamilto- 
nian J2]) and, therefore, can be considered as a dynamical 
counterpart of the latter. Like the QHF order parameter 

(|16p . this mass term completely breaks the Z^y symme- 
try. 

It is important that in both monolayer and bilayer 
graphene, these two sets of the order parameters neces- 
sarily coexis t 26 ' 27 ' 37 and are produced even at the weak- 
est repulsive interactions between electrons (magnetic 
catalysi o 10-12 ). The essence of this phenomenon is an ef- 
fective reduction by two units of the spatial dimension in 
the electron-hole pairing in the LLL with energy E = 0. 

Beside these parameters, the Green's function includes 
also the chemical potential ^ related to the charge den- 
sity S{ s (^\ s^C.s)- Because of that, it will be conve- 
nient to consider their combinations m = fx + fJ-3 and 
fil = /1 — (13. Let us emphasize that all the parameters 
(fi s , pL s , A s , and A s ) should be viewed as operators with 
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well defined values only when projected onto quasiparti- 
cle quantum states, in which case they become functions 
of the Landau level index n: fi ns , fins, A ns , and A„ s . The 
derivations and explicit expressions of the full / free quasi- 
particle Green's functions are presented in Appendix [X] 
Concerning the full Green's function, it is interesting 
to note that the order parameters in the lowest Landau 
level enter the quasiparticlc energies only through two 
independent combinations^ Ef na = — (/i„ s + A ns ) + 



£,{fas — A„ s ), where n = 0, 1. This is the consequence of 
the spinor structure in LLL wave functions, which have 
only half of the components nonvanishing (the valley and 
layer degrees of freedom are not independent in the LLL). 
Formally, this is also the reason why the roles of QHF and 
MC order parameters are indistinguishable in the LLL 
approximation. Without the loss of generality, we will 
assume in the following that only A„ s and A„ s are non- 
trivial parameters in the lowest Landau level (n = 0, 1). 



According to Eqs. (|A10[) and (|A16[) in Appendix [3] the Green's function and its inverse contain the same overall 
Schwinger phase, which describes the breakdown of the conventional translation invariance in a magnetic field. After 
substituting these functions into the gap equation ([5]), omitting the phase on both sides, and performing a Fourier 
transform with respect to the time variable, we arrive at the following equation for the translationally invariant part 
of the Green's function: 



2?r / 2tt 



= S-V;r)- / ^ / ^J ( 

1 f dto' 

2 J 2^ 
i f dui 1 

27 



pdp 
~2T 



Jo{pr) G(Ly ; r)-T 3 G(cy;r)T 3 - «;p) 



[Pi 



(0 - V 2 (0) tr {[Pi(0 - P 2 (0] G(u'; 0)} V^ c (0)S(r) . 



(19) 



After substituting the explicit form of the translationally invariant functions given in Eqs. (|A11[) and (|A17I) . we 
calculate the traces and obtain equations for the order parameters n, p,, A, and A. 

Furthermore, assuming these order parameters are energy independent, we can analytically integrate over the 
energy and momentum on the right hand side of the gap equation (|19p . Finally, after projecting out onto different 
orbitals (associated with different Lagucrre polynomials), we arrive at a final set of algebraic equations for the order 
parameters in all Landau levels: 



E £0s — Ms - £^0 



hide 



-I ofi (E^ 0s ) S ign(E^ s ) - h, (E^ s ) sign(££ s )] 



TlUJ c ^ — ^ . \ / ( ^6n f s 

In' fl(Mtn> s + Pen's) sign(M^„/ s + fan's) 1 + T7 

n'=2 V M «"' S 

E tn'AMtn's - fan's) sign(M in > s - fan's) fl 
n'=2 



M 6r 



1 1 ™ max A 

-sign(^ 0s ,) + -sign(^ ls ,) - W^f: 6 ( M k's' - 4n's') 



(20) 



Tiuj 



TlUJr 



TlJjJr 



) sign(M 5 „- s + fan's) 1 



n'=2 



t±£n's 
M in 's 

A^„' S 

M in , s 

-sign(4 0s ,) + -sign(4 ls ,) - ]T J^ml'n's' ~ 4n's'] 

n 1 —0, ^ 



) sign(A^t n / s - fan's) 1 



(21) 
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for n — and n = 1, respectively, and a pair of equations, 



^ ns - fi s + A ina - £A = _ £ [-/ 0in (££ J s ign(££ 0s ) - I liB (££ a )sign(££ l5 )] 



) sign(M^„/ s + /if n ' s ) ( 1 + 



n'=2 



~T E ^'^(M^n's - H£n's) sign(M ?)l / s - H£n's) ( 1 - 



n'=2 

TliUJ 



A^n's 



E* 



1 1 max A 

±sign(if, 0s ,) + isign(^ ls ,) - £ (M|„ 

n—2 s 



2 \ 



(22) 



/Ugns ~ (J-s ~ A^ ns + £A = — £ ^ / n '-2,«-2(M ? „' s + /if„' s ) sign(M ? „< s + fi^n's) ( 1 



n'=2 



TlLO r 



^ / n '-2,n-2(A% l ' S - |"£n' S ) sign(M ?n < s - //£„< s ) ( 1 + — 



A^n's 

At-/, 



E? 



1 • '-iL \ , 1 • / t^L \ Af/ n / s / 

1 rfflgn^!,,) - 2^ 



-sign(£|; 0s 



n'=2 



6(M? 



2 \ 
Vgn's'l 



, (23) 



for each n > 2. Here, p, s = /i — sZ is the effective chemical potential that includes the shift due to the Zeeman 
energy (s = ±1 corresponds to up and down spins, respectively), fi is the chemical potential itself, and «il = 
e 2r )\d/ (2h 2 v 2 7 n) tn 0.354/k is a dimensionlcss interlaycr coupling constant. 



The other notations are 



the following functional dependence on the energy: 



ns • 



(24) 
(25) 



\ns = ^A\ ns + (?iw c ) 2 n(n - 1), for n > 2,(26) 



E, 



£ 1 1 s 



-A ns -£A ns , for n = 0, 1. 



(27) 



It is shown in Appendix |A] that while i?L s in Eq. ([27]) 
are the quasiparticle energies in the two lowest Landau 
levels, the quasiparticle energies in higher Landau levels 
(n > 2) are given by the following expression: 



% s = -M£n«±M£ ns , for n > 2. 



(28) 



The set of gap equations (j2T)|) - (f2"3"]) determines the 
order parameters /it, /2, A, and A as functions of the 
Landau level index n and the external parameters: the 
chemical potential po, the magnetic field B, and the ap- 
plied electric field E±, which is expressed through the 
top-bottom gate imbalance Ao, E± = 2Aq/ (ed)X The 
energy-dependent coefficient functions I n t^ n (E) in these 
equations are given in Eq. (|B1[) in Appendix [B] Each of 
the dimensionless coefficient functions I n '.n(E) can be 
calculated numerically with a relative ease. However, 
when repeatedly solving a complete set of gap equations 
with many Landau levels, it is greatly beneficial to use 
cither a tabulated set of these functions or approximate 
analytical expressions. It appears that the numerical re- 
sults for functions I n i^ n (E) can be well approximated by 



In',n(E) 



E 



1 + d n ',n hul 



(29) 



which is usually valid in a wide range of energies (0 < 
E < 40Huj c ). In the case of B = 2T and K = 2, for 
example, the corresponding coefficients for < n',n < 10 
are summarized in Table IIIII in Appendix [B] We show 
only the upper triangular part of the corresponding tables 
because all coefficients arc symmetric with respect to the 
exchange n' and n. 

The approximation with static screening can be eas- 
ily obtained from the more general dynamical one in 
Eq. (|B1[) by substituting E = and taking hi = 
(i = 1,4). The corresponding numerical calculation be- 
comes very easy and one can tabulate the static coeffi- 
cients I n ',n on the fly. 

Before concluding this section, let us also briefly dis- 
cuss the truncation of a formally infinite set of the gap 
equations. From a physics viewpoint, it is clear that the 
low-energy effective model is valid only in a finite range 
of energies (up to about energy A « 71 /4 ~ 1000 K)J- 
Therefore, the sum over Landau levels in Eqs. (|20|) - ([23]) 
should be truncated at n max ~ fcsA/(?iu; c ) ~ 40/ (B [T]), 
i.e., n max is different for different values of B and it de- 
creases with increasing B. The use of such a prescription 
implies that the model at hand is appropriate only for 
magnetic fields B < 40 T. In order to consider stronger 
magnetic fields, B > 40 T, or include Landau levels be- 
yond n max , one should utilize the microscopic four-band 
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niodclii This is beyond the scope of the present work. 
In the opposite case of weak magnetic fields, B < 1 T, 
the model becomes expensive numerically and should be 
treated differently. In particular, the limit B — > (no 
magnetic field) can be described only if an infinite num- 
ber of Landau levels are taken into account, which would 
require a different approach to the problem. 

IV. NUMERICAL RESULTS 

In this study we concentrate on the ground states at 
the neutrality point. For the purposes of our analysis, 
it is sufficient to fix the value of the chemical potential 
jUo = 0. We will see that there are two main solutions to 
the gap equation: the spin- and layer-polarized ones4£ In 
order to determine which of them corresponds to a true 
ground state of bilayer at given values of gates bias and 
magnetic field, we must compare their free energy den- 
sities. The corresponding expression for the free energy 
density is derived in Appendix [Cl 

A. Landau level mixing in static screening 
approximation 

Before we proceed to a more complicated analysis 
of the gap equations with dynamical screening in Sub- 
sec. IIVB1 let us first present our main results for the 
static screening case. Such an approximation is instruc- 
tive for getting a better understanding of the Landau 
level mixing effects. 

Let us start our analysis by considering a simple bench- 
mark case with B = 2T (n max = 20) and k = 2. For a 
fixed value of the voltage imbalance Ao , we typically find 
two types of competing solutions: spin-polarized (SP) 
and layer-polarized (LP) ones. The examples of both 
types of solutions are shown in Fig. [5J which are ob- 
tained at a fixed value of Ao ~ 0.3?ia; c , which is close to 
the critical value Ao, C r ~ 0.32tkj c separating the SP and 
LP phases (sec Fig. [3] and its discussion below). Note 
the different energy scales on the left and right panels in 
Fig.H 

As one can see from this figure, the QHF parameters 
(fx ns , fins) and the MC ones (A„ s , A„ )S ) coexist in all 
Landau levels. The Haldane mass A ns and the chemical 
potential fj, ns have different signs for up and down spins 
in the SP solution. The values of |Ao, s | — 0.83fioj c and 
| Ai a | — 0.74?iw c in the LLL are essentially larger than 
those for higher LLs with n > 2, which arc slowly de- 
creasing with increasing n: from |A2 lS | — 0.25?k*; c down 
to |A 2 o, s | — 0.137iw c . Its QHF counterpart fx ns behaves 
similarly: \fi2,s\ — 0.147kJ c and |/X20,s| — 0.11fik) c . (Recall 
that hu c ~ 430 meV for B = 2T.) 

As to the dynamical voltage imbalance (Dirac mass) 
A ns , it is suppressed with respect to the bare voltage 
Ao ss 0.3hu> c used in this figure: Ao. s ~ 0.16huj c and 
Ai iS ~ O.I6?ia; c in the LLL, while A2 iS ~ 0.23?k*; c and 



A20.S — 0.l5tkj c . The value of its QHF counterpart fl ns is 
very small. It starts from p,2,s — 0.009?icl> c and decreases 
down to the values of order 10~ 4 huj c at large n. Thus, 
taking into account the dispersion relations in Eqs. (|27p 
and ([28]) . we conclude that, as expected, the splitting of 
the levels with opposite spins is responsible for generating 
a gap in the SP solution (see Fig. [3] and its discussion 
below) . 

In the LP solution, the values |A ns | and |/x ns | are small. 
In fact, while |Ao s | and |Ai s | are equal to the Zeeman 
energy Z = Q.027Tiu> c , all A ns with n > 2 vanish. The 
chemical potential \[i ns \ is equal to the Zeeman energy 
for all n. As to the parameters A ns and jl ns , the val- 
ues of the former in the LLL are significantly larger than 
those in the higher LLs: Ao, s = 0.68Tioj c , Ai jS = 0.59?iu; c , 
while A 2>s = 0.07?iu; c and A 20 , s = 0.007^o; c . Its QHF 
counterpart parameters /x n>s are /t 2 ,s = —0.12huj c , = 
— 0.10?k*; c , and /U2o,s = — 0.08?ia; c , whose absolute values 
are considerably larger than the Zeeman energy. From 
the dispersion relations in Eqs. (|27|) and (|28p. we con- 
clude that the splitting of the levels assigned to different 
valleys, K (£ = 1) and K 1 (£ = —1), is responsible for 
generating a gap in the LP solution (recall that the valley 
and layer degrees of freedom are not independent in the 
LLL). 

Although the values of the QHF and MC parameters 
are different at other values of the bare voltage imbal- 
ance Ao, the main characteristics of their dependence on 
the Landau index n remain qualitatively similar. Instead 
of showing the QHF and MC parameters as functions of 
Landau level index n for other values of Ao, it is conve- 
nient to summarize the Ao dependence of the SP- and 
LP-type solutions by presenting the spectra of the first 
few low-energy states in Fig. [3j The first two panels show 
the energies of the first few Landau levels for the SP and 
LP solutions. The free energies are compared in the right 
panel of the same figure. 

As we see, the spectrum of the LLL with n = and n = 
1 is qualitatively different from that of the n = 2 LL. The 
roots of this difference are in the form of the spectrum in 
the bilayer model without interactions^ While the ener- 
gies of LLL states E^ 0s and E^i s equal zero, there are pos- 
itive and negative energy bands E^ ns = ±y/n(n — T)hu> c 
for each of the higher LLs with n> 2. 

The SP solution is the ground state at small values 
of Ao, while the LP one becomes the ground state at 
large values of Ao. The first order phase transitio n 27 ' 31 
from one regime to the other occurs at the critical value 
Ao, C r ~ Q.32hu> c . In terms of the applied electric field, 
this is equivalent to E± tCT w 7.87mV/nm (at fixed value 
B = 2 T), which is somewhat smaller than typical values 
measured in the experiment . 15 ' 17 " — 

It is instructive to compare the gaps in the energy spec- 
tra for both the SP and LP solutions at a fixed value of 
magnetic field, B = 2 T. At Ao = the energy gap for 
the SP solution is E^ p = El P 1A - E^^ ~ 1.49hu c , 
which is considerably larger than the corresponding gap 
for the LP solution, E^ p = £^ p 14 - £^p t ^ o.80?ia; c . 
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FIG. 2: (Color online) The dependence of the dynamical parameters on the Landau level index n for the SP (top) and LP 
(bottom) solutions at a fixed value of Ao ~ 0.3huj c in the approximation with static screening. The magnetic field is B — 2 T. 



At the critical value Ao. cr ~ 0.32hus c , however, both 
gaps happen to be approximately equal to l.lAhuj c . 
These values are also comparable to those in some 
experiments ) 15 ' 16 but are considerably smaller than the 
recent experimental values in suspended bilayer graphene 
with a high mobility.— 

From a theoretical viewpoint, it is interesting to in- 
vestigate the dependence of the main results on n max , 
which plays the role of the high-energy cutoff parameter. 
While the qualitative competition between the SP and 
LP states remains the same, the critical values E± tCT de- 
crease slightly with decreasing the number of Landau lev- 
els (i.e., shrinking the available phase space). The actual 
critical values of the applied electric field (at B = 2 T) 
for several different values of the Landau level cutoff pa- 
rameter n max are given in Table HI These data show that 
the critical value of -Ej_ jCr may have a sizable uncertainty, 
associated with the choice of n max . Indeed, its approx- 
imate physical value is determined only semi-rigorously 
as the number of Landau levels below the energy cutoff 
A~7i/4. 

Now, let us discuss the dependence of our numerical 
results on the magnetic field. Neglecting a weak depen- 
dence of the interaction coefficients I n ,n' on the mag- 
netic field B in the static approximation, one may claim 
that the full dependence of the dynamical parameters on 
B can be completely restored from simple scaling argu- 
ments. As seen from the gap equations (J2DJ) - (J2SJ), the 
main dependence on the magnetic field comes through 



the overall factor huj c on the right hand side of all equa- 
tions. (Recall that, in the static approximation, the co- 
efficient functions I n ',n are energy independent.) This 
dependence on hu c can be removed by introducing di- 
mensionless energy parameters, measured in units of hu> c . 
If this were the only dependence on the magnetic field, 
the results for all dimensionless ratios, such as /i„ s /?iw c , 
fins/fiWc, A ns /huj c , A ns /huj c and all energies E ins /huj c 
would be the same for any value of B. 

An additional weak dependence on the magnetic field 
enters the gap equations indirectly through the cutoff pa- 
rameter n max and through the x-parametcr, defined after 
Eq. ([8]). The latter appears in the definition of the coef- 
ficient functions I n ', n even in the static approximation. 
However, the corresponding dependence on the magnetic 
field is relatively weak. When the value of B increases by 
two orders of magnitude (from 0.1 T to 10 T), the numer- 
ical values of all off-diagonal coefficients (n' ^ n) 
change only by a few percent. There is a stronger de- 
pendence on the field in the diagonal coefficient I n . n , but 
even this one changes only by about a factor of 2. 

In order to see the deviations from the linear scaling 
laws, we did perform a careful numerical dependence of 
the results on the magnetic field. Note that, in such a 
study, one has to adjust properly the maximum value of 
the Landau levels kept in the gap equations. Recalling 
that n max = fcsA/(7kj c ) ~ 40/(5 [T]), we will consider 
several values of magnetic fields (e.g., 2T, 4T, 5T, 8T) 
which lead to integer values of n max . The best fits to the 
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FIG. 3: (Color online) Quasiparticle energies in the lowest three Landau levels, n = (solid lines), n = 1 (dashed lines) and 
n = 2 (dotted lines), for the SP (left panel) and LP (middle panel) solutions obtained in the approximation with static screening 
at B — 2T. Colors of the lines correspond to specific values of quantum numbers (£, s): red to ( — ,1), green to (— , f), blue to 
and purple to (+,t)- The right panel shows the free energies of the same solutions as functions of Aq. 



TABLE I: Critical values of Ao and E± for several choices of the cutoff parameter n max in the approximation with static 
screening at B = 2 T and « = 2. 
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Ao,cr /huJc 


0.204 


0.234 


0.250 


0.260 


0.268 


0.275 


0.293 


0.318 


-Ex.cr [mV/nm] 


5.02 


5.75 


6.14 


6.40 


6.60 


6.75 


7.20 


7.83 



critical lines are given by the following expressions: 

K = 1 : Ejl, C t ^ (3.82 + 5.985 [T]) mV/nm, (30) 
K = 2 : E x , a ~ (1.69 + 3.075 [T]) mV/nm. (31) 

One of the qualitative deviations from the linear scal- 
ing is the appearance of a small positive offset in the 
dependence of the critical value E±_ cr on the magnetic 
field. It is noticeable that such an offset would be absent 
without Landau level mixing. Although strictly speaking 
we cannot consider the limit of vanishingly small mag- 
netic fields using the present approach, the interpolation 
of this result qualitatively agrees with the experimental 
dat a 15 ' 17 " — showing that the SP state extends all the 
way to B = at sufficiently small electric fields. This 
seems to suggest that it can be the ground state of bilayer 
graphene at the neutral point in the absence of external 
fields. 

Comparing these results with the earlier studies with- 
out Landau level mixing, we see that the size of the 
energy gaps in the SP and LP type ground state as 
well as the critical values of the applied electric field 
substantially increased. However, even the inclusion of 
Landau mixing did not resolve completely the discrep- 
ancy with the corresponding experimental values in clean 
samples^ which appear to be still larger than our pre- 
dictions in a model with static screening. It is critical, 
therefore, to go beyond the static approximation. 

B. Dynamical screening 

Now let us turn to the analysis of the gap equations 
with dynamical screening. Computationally, this is much 



harder than the static case. It is a nontrivial energy de- 
pendence of the coefficient functions I n i ^ n (E) that makes 
the numerical computations considerably slower. Also, 
a highly nonlinear nature of the gap equations makes 
finding numerical solutions much harder. Before even 
solving the gap equations, in fact, one first needs to com- 
pile all coefficient functions I n ^ n (E) with n',n < n max 
using their definition in Eq. (|B1[) . Fortunately, the task 
is somewhat simplified by the observation that each of 
these functions can be fitted quite well with a simple 
Pade approximant of order [2/1], see Eq. (|2"9")l . The nu- 
merical values of the coefficients for the fits of I n > , n {E) 
with n' ,n < 10 are listed in Table IILTI in Appendix [Bl 

Let us start from considering some key results in the 
same benchmark case of a fixed magnetic field, B = 2 T. 
By reviewing the functional dependence of the coeffi- 
cients I n i^ n (E), compiled numerically, we find that these 
functions grow substantially with energy. For example, 
typical values of the coefficients I n ',n(E) with n' , n < 10 
at E = huj c appear to be about 20% to 90% larger then 
their values at E = (i.e., static screening approxima- 
tion). From this fact alone, it is natural to expect that 
the dynamical parameters as well as the associated en- 
ergy gaps in the spectra should be larger in the approx- 
imation with dynamic screening. This is indeed what 
we see. The examples of the SP and LP solutions are 
shown in Fig. [4j which are obtained at a fixed value of 
Ao ~ 1.02?ia; c ~ 4.39 meV that is close to the critical 
value Ao jC r — 0.99?kj c (note the different energy scales 
on the left and right panels in Fig. 2]). This figure illus- 
trates that, as in the case of the static screening, the MC 
and QHF dynamical parameters necessarily coexist. 

Although the MC and QHF parameters in the case of 
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FIG. 4: (Color online) The dependence of the dynamical parameters on the Landau level index n for the SP (top panels) and 
LP (bottom panels) solutions at a fixed value of Ao « hu c and re = 2. The numerical results are obtained in the approximation 
with n max = 20 and dynamic screening. The magnetic field is B — 2T. 



the dynamical screening are considerably larger than in 
the case of the static screening, their functional depen- 
dence on the Landau index n and some other important 
characteristics are similar. As we see from Fig. HI the 
Haldane mass A„ s and the chemical potential /i ns have 
different signs for up and down spins in the SP state. 
The values of |A , S | ~ 2.20fiw c and |A 1;S | ~ 2.02^w c 
in the LLL are substantially larger than those in higher 
LLs with n > 2, which are slowly decreasing with n 
from |A2, S | — l-20?ia; c down to |A2o. s | — 0.40?ia; c . Its 
QHF counterpart [i ns behaves similarly, decreasing from 
\fJ,2,s\ — 0.35hui c to |/X2o,s| — 0.15fi.w c . 

As to the dynamical voltage imbalance (Dirac mass) 
A ns , unlike the static case, its values in the LLL and 
n = 2 LL are not suppressed with respect to the bare 
voltage Ao « 1.02?kJ c . Namely, Ao, s — 1.02hu> c and 
Ai iS ~ 1.07Tiuj c in the LLL, while A2 iS ~ 1A3Juj c . Its 
value at the largest n = 20 is A 2 o, s — Q.54Tioj c . As in 
the static case, the values of its QHF counterpart jl ns 
are much smaller, starting from //2, s — 0.lOhui c and de- 
creasing down to the values of order 10 _3 7ia; c at large 
n. Therefore, we conclude that the splitting of the levels 
with opposite spins is responsible for generating a gap in 
the SP solution (see also Fig. [5] and its discussion below). 

In the LP solution, the values |A ns | and \jJL ns \ are 
small, although unlike the static screening case, their 
values in the LLL are larger that the Zeeman energy 
Z = 0.027?^. In fact, |A , S | ~ 0.06So; c and |A M | ~ 



O.Q6Tiuj c , while all A ns with n > 2 are of the order 
10 _3 ?ia; c . The chemical potential \/J, ns \ slowly decreases 
from \n2,s\ — O.Q55huj c to \fJ.20,s\ — 0.001?ia; c . As to 
the parameters A ns and p, ns , the values of the former 
in the LLL are significantly larger than those in the 
higher LLs: Ao, s — 2.08?ia; c , Ai. s ~ 1.93?ia; c , while 
A2, s ~ 1.36?ia; c and A20.S — 0-30fioj c . Its QHF coun- 
terpart /2 n s is pL2. s — — 0.25?iuj c , /i3 jS ~ —0.2Atkj c . and 
p,2o,s — —0-12huj c , whose absolute values are consider- 
ably larger than the Zeeman energy. From the dispersion 
relations in Eqs. (|2"T|) and (|2"8|) , we conclude that the split- 
ting of the levels assigned to different valleys, £ = ±1, is 
responsible for generating a gap in the LP solution. 

The energy spectra of the low-energy states are pre- 
sented in Fig. [5] A few comments are in order here. In 
addition to the expected large values of the gaps, we also 
find that (i) the energies have a rather complicated de- 
pendence on the voltage imbalance Ao, which substan- 
tially deviates from a linear dependence, (ii) there are 
several points of level crossings and (iii) a nonstandard 
order of the Landau levels with the lowest energy state 
being the n = 2 Landau level. The reason of the last phe- 
nomenon is connected with the following feature: As one 
can sec in Fig. 01 the values of some dynamical parame- 
ters in the LLL are much larger than those in the n = 2 
Landau level. Interestingly, this feature takes place for 
both the SP and LP solutions. 

The energy gap in the SP ground state is given by 
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FIG. 5: (Color online) Quasiparticle energies in the lowest three Landau levels, n = (solid lines), n = 1 (dashed lines) and 
n = 2 (dotted lines), for the SP (left panel) and LP (middle panel) solutions obtained in the approximation with dynamic 
screening at B = 2 T. Colors of the lines correspond to specific values of quantum numbers (£, s): red to ( — , l), green to ( — , t)i 
blue to (+,4), and purple to ( + ,t)- The right panel shows the free energies of the same solutions as functions of A . 




E g*p = E -a-X ~ E +a,-\- At A = 0, the corresponding 
value of the gap is ~ AA2hw c at B = 2 T. This gap 
grows with increasing magnetic field, and the correspond- 
ing dependence is approximately a linear function of B. 
Our numerical results and linear fits are shown in Fig. [6] 
for the cases of dynamical as well as static screening. 

At the critical point, Ao, cr ~ 0.997kj c , the SP gap E^ p 
and the LP one, Ej£L = -E+,2.f — -E-,2,t> take different 
values: E^?„ ~ 1.95huj c and E^ ~ 3.24hw c . In other 

gap gap 

words there is a jump in the gap at the critical point 
shown in Fig. [7] For comparison, the gap in the case 
of the static screening is also shown in this figure. As 
one can see, at the critical point it has a kink that is 
smoother than a jump singularity. As discussed in Sec.lVl 
below, the presence of such singularities at the critical 
point is relevant for understanding the behavior of the 
conductivity observed in experiment^ 

Another consequence of the approximation with the 
dynamical screening is a substantial enhancement of the 
critical value of the applied electric field, at which the 
phase transition from SP to LP states occurs. For the 
results presented in Fig. [SI for example, the critical value 
is -B^.cr ~ 24.31 mV/nm (as seen from the figure for 



the free energy, the corresponding voltage imbalance is 
A ,cr ~ 0.99Tllj c ). In order to appreciate the effect of 
Landau level mixing, we present the numerical values of 
Ao and E± for several choices of the cutoff parameter 
n max in Table ILT1 

The corresponding numerical results for magnetic field 
dependence of the critical value of E± are shown in Fig. [8] 
Just like in the static case, we adjust the maximum value 
of the Landau levels as follows: n max = kB-h/huj c (B), 
where A = 1000 K is a fixed cutoff. The data can be well 
approximated by the following linear dependence with a 
nonzero offset: 



E± cr ~ (19.85 + 2.135 [T]) mV/nm, (A= 1000K). 

(32) 

Comparing this expression with that in Eq. (|31[) , we con- 
clude that the dynamical E±, cr is considerably larger 
than the static one. 
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TABLE II: Critical values of Ao and E± for several choices of the cutoff parameter n max in the approximation with dynamical 
screening at B = 2 T and k — 2. 





1 


2 


3 


4 


5 


6 


10 


20 




0.241 


0.336 


0.412 


0.476 


0.530 


0.578 


0.732 


0.988 


E x ,ci [mV/nm] 


5.92 


8.27 


10.2 


11.7 


13.0 


14.2 


18.0 


24.3 




2 4 6 8 

BIT] 



FIG. 8: (Color online) The dependence of the critical value 
of E± on the magnetic field in the case of « = 2. 



V. DISCUSSION 

The two ingredients, the Landau level mixing and the 
dynamical screening, play an important role in the dy- 
namics in bilaycr graphene in a magnetic field4£ In fact, 
their role is more important than that in monolayer 
graphene (compare with Ref. l42h . This is a reflection 
of the general fact that the Coulomb interaction in bi- 
layer graphene plays a more profound role. The latter is 
the consequence of a much smaller characteristic energy 
scale huj c ~ 2.15B [T] meV in bilayer graphene as com- 
pared to the Landau energy scale si ~ 26 -y/ ' B [T] meV in 
monolayer. 

In this study, we derived and solved a complete set 
of gap equations with Landau mixing in the Hartree- 
Fock approximation with the frequency-dependent polar- 
ization function calculated at the charge neutrality point. 
The competition between the spin-polarized and layer po- 
larized states is studied in detail by varying the applied 
electric field (or, equivalently, the voltage imbalance be- 
tween the top and bottom layers) and the strength of the 
magnetic field. 

It was found that the critical value of the applied elec- 
tric field is a linear function of the magnetic field, i.e., 
E±,a = + aB, where is an offset electric field 
and a is the slope. The offset electric field and energy 
gaps substantially increase with the inclusion of dynam- 
ical screening compared to the case of static screening. 
The solutions to the gap equations clearly demonstrate 
that the QHF and MS order parameters necessarily co- 
exist. 

We demonstrated that the use of the static screening 



approximation badly underestimates the strength of the 
direct Coulomb interaction. Taking into account the ef- 
fects of dynamic screening leads to energy gaps that are 
about a factor of 2 to 3 larger than those in the approxi- 
mation with static screening. The increased values of the 
gaps as well as more pronounced nonlinearity of the gap 
equations result in a rather non-trivial energy spectrum. 
In particular, we find that the energy of the n = 2 Lan- 
dau level becomes smaller than the n — 0, 1 ("lowest") 
Landau levels. 

Let us now compare the theoretical model results with 
the experimental data in Refs. Il5lll6lfl9l . One of the most 
noticeable features of the present model is the approxi- 
mate linear dependences of the energy gap E gap and the 
critical electric field E±, CI on the magnetic field. We find 
not only that these functions are linear in B, as many 
previous model calculations showed, but also that they 
have a nonzero intercept (offset), which appears as a re- 
sult of the Landau level mixing. 

The results for the SP gap (at Ao = 0) as a func- 
tion of the magnetic field are shown in the left panel in 
Fig. [5] As we see, the results in the case of static screen- 
ing are in a reasonable agreement with the experiment 
in the low-mobility graphene samples^ The results for 
the gaps in the case of dynamical screening are in or- 
der of magnitude agreement with the experiment in the 
high-mobility graphene samples*^ On the other hand, 
the slopes and the intercepts of theoretical lines are not in 
excellent quantitative agreement. There may exist many 
potential reasons for the discrepancies. One of them is 
the limitations of the model itself, for example, the use 
of the zero width of the Landau levels. The other is the 
approximations used in the analysis of the model. Per- 
haps the most important limitation of the second type is 
the energy independent ansatz for the gap parameters. 

As to the dependence of the critical electric field E±, CI 
on the magnetic field in Eq. (f3"2"|) , we find that it is a linear 
function of B just like in the experimen t 15 ' 19 and that the 
value of the offset is in a reasonable agreement with the 
experimental values (see the right panel in Fig.[S|). 

The theoretical slope 2.13mVnm~ 1 T~ , however, ap- 
pears to be much smaller than the experiment value, 
12.7mVnm _1 T — . This discrepancy may have its roots 
in disorder, which presumably plays a profound role in 
real samples. One may suggest, for example, that an 
external electric field is more effective in clean samples 
(considered in our model) and, therefore, its critical value 
E^. cr should be smaller than that in samples with impu- 
rities. This conjecture agrees with that E±^ cr is smaller in 
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FIG. 9: (Color online) Comparison of the energy gap -E gap in the SP state as a function of the magnetic field at E± = with 
the experimental results of Refs. [Ir3ll9l (left panel). Comparison of the dependence of the critical value of E± on the magnetic 
field with the experimental results of Refs. Il5lfl9l (right panel). 



the experiment^ with high-mobility samples than that 
in the experiment^ with low-mobility ones (see the right 
panel in Fig. |9] above). 

It is appropriate to mention another key experimental 
observation in Ref. [ID that the two-terminal conductance 
is quantized except at particular values of the electric 
field E±. For the v = state, the experimental data 
shows that the conductance is quantized except for two 
values of E± Let us argue that these two values corre- 
spond to E±_ = ±_E_l jCI ., where E±^ cr is the critical value 
at which the phase transition between the SP and LP 
phases takes place (see Fig. where we consider only 
nonnegative values of Ao = eE±d/2). As one can see 
in this figure, the gap A„ = o has a maximum at E± = 
and a minimum at E± = E± cr . This implies that while 
the conductivity is suppressed almost everywhere, it is 
enhanced for the two values of E± = ±-Ej_ jCr (compare 
with Fig. 2C in Ref.[lj|. 

In this paper, we did not investigate the limit of weak 
magnetic fields, which is outside the scope of the present 
paper. In fact, the numerical approach used here is not 
efficient in such a limit because the number of Landau 
levels becomes too large and the importance of level dis- 
cretization itself diminishes. It is still interesting to spec- 
ulate about the physical meaning of our result when ex- 
trapolated to vanishing fields. A non-zero value of the off- 
set electric field E'f 1 , needed in the fit of the critical line, 
suggests that the ground state of bilayer graphene with- 
out magnetic and electric fields can be a spin-polarized 
(ferromagnetic) state. In the future studies, it will be in- 
teresting to address this question in detail by considering 
the weak field limit. 

In the recent experimental papcr^ a phase transition 
to the nematic state was observed in bilayer graphene 
without magnetic field. The nematic state has been pre- 
dicted in theoretical works^ It breaks the rotational 
symmetry and keeps quasiparticles gapless.— It would be 
interesting to study a competition between the gapped 
state and the nematic state due to the Coulomb inter- 
action and other interactions in bilayer graphene in a 



magnetic field. 

In the future, it would be also interesting to extend 
our study to the case of a more general model of bilayer 
graphene formulated in terms of the four-component 
spinors, which has a much larger energy range of ap- 
plicability and allows to address the competition of the 
SP and LP states in the very strong magnetic fields. The 
dynamics at larger energies should reveal an interesting 
cross-over to the regime when the top and bottom layers 
of bilayer effectively decouple. 

It will be also interesting to extend the analysis of this 
paper to the case of quantum Hall states with nonzero fill- 
ing factors. It is natural to expect that the Landau level 
mixing and dynamical screening will again substantially 
modify the corresponding results obtained in the LLL 
approximation . 2 1 ' 25 ~ £L Finally, there remains the ques- 
tion whether there are nonuniform states, e.g., such as 
the helical and electron crystal states^ 2 - that are ener- 
getically more favorable in bilayer graphene under certain 
conditions. 
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Appendix A: Quasiparticle Green's function 

The quasiparticle Green's function can be obtained in a standard way by making use of a complete set of eigenstates 
for the low-energy Hamiltonian in bilayer graphene in a magnetic field. Schematically, 



Sr( W ; P> r') = X; 

{ki} 



w - E {ki} 



(Al) 



where {ki} stands for a complete set of quantum numbers that uniquely define the eigenstates. Here and in the rest 
of appendixes we put Ti = 1. 

In the Landau gauge, A = (A x , A y , A z ) = (0, B±x, 0), the corresponding set of eigenfunctions for the free Hamilto- 
nian reads 



*fc,£,n,T=-l(r) 







n = 0,l, 



**» JM n + T^A u n - 2 {x) 




V^Wn V T jM n -T^A U n ( X ) 



n > 2, 



(A2) 
(A3) 



where x = fc£ + M n — y A 2 , + uj 2 n(n — 1), r = ±, and 

1 e~x 2 /2 

«n(x) 



(A4) 



The corresponding energy eigenvalues are ±M„. 

By making use of the above complete set of eigenstates, it is straightforward to derive the following explicit form 
of the free quasiparticle Green's function: 



S(u; r, r 



S(uj;r-r') 



n max 

+ E 



E 



(w + /i 8 - ^A )L„(2;) + +^A )L„- 2 (2)„ 

■ r - H ; — tt; ~ To r 4 



n=2 

_ / r 2 _ 
2[(u; + p s f-M^ \rl 



Wat*)!, 



(A5) 



(A6) 



where 2 = r 2 /(2£ 2 ), fj, s = fio — sZ, r± = x ± and "P± = (1 ± ts)/2. Note that the translational invariance of this 
function is spoiled only by the Schwingcr phase factor e**( r ' r ). I n the Landau gauge used here, the explicit form of 
the phase reads 



(x + x')(v — v') , „ n 

Using the same approach, we can also derive a similar representation for the inverse Green's function: 
S-Vir.rO = e i *<- r > r ' ) S- l {w;r-r'), 



S _1 (w;r) 



2tt£ 2 

n max 

E 

n=2 



(w + /*, + fA )[£o(z) + Li(«)] P_ 



(w + jl, + £A )L„(z)7>_ + (w + /2 S - £A )L„_ 2 (z)7> + - ^ ( ° 2 ^ j L 2 _ 2 (z) 



(A7) 



(A8) 



.(A9) 



Let us emphasize that the Green's function and its inverse have the same Schwinger phases. One can also show that 
the phase remains exactly the same even for the full quasiparticle Green's function (see below). This property plays 
an important role in the derivation of the gap (Schwingcr-Dyson) equation, which takes a very simple form after the 
same overall non-zero Schwinger phase on both sides of the equation is eliminated. 
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The full Green's function, which incorporates the effects due to dynamically generated order parameters of quantum 
Hall ferromagnetism and magnetic catalysis, can be derived in the same way as the free function. The final result 
reads 



G(ui; r, r 



:/2 



G(w;r-r'), 
Lo(z) 



2ir£ 

n max 

E 



V- 



r 2 



^2[{w + nZnsY-Ml na ]P \r\ 



(w + ^, ls -A^ s )i n (z) (uj + n ins + A^ ns )L n ^ 2 {z) 

■ r- H : : rr; n 



(w + ^ ns ) 2 -M, 2 



(w + // e „ s ) 2 - M, 2 



where 



£A ns , for n = 0, 1, 



M c „ s = ^A| n8 +w2n(n- 1), for n > 2. 

The corresponding inverse Green's function is 

G-\uj;r,r') = e^G-^w; r - r'), 
e -»/a r 



G- x (o;;r) 



(cj - E^ a )L (z)V- + E^ s )L x {z)V- -Y,ik 



r 2 



2£ 2 I r 2 



^ [(w + Mens + A ?ns )L n (2)"P_ + (w + fj,^ na - A ins )L n ^ 2 (z)'P+} 



Ll- 2 {z) 



(A10) 



(All) 



(A12) 
(A13) 
(A14) 

(A15) 



(A16) 



(A17) 



While E^ ns are the quasiparticle energies in the two lowest Landau levels, the quasiparticle energies in higher Landau 
levels (n > 2) are given by the following expression: 



E, 



-fi^ns ± M£ na , for n > 2. 



(A18) 



It should be emphasized that, while the Green's functions G(u; r, r') and G r, r') are inverse of each other, their 
translationally invariant parts, G(ui;r) and (w;r), arc not. Another important property, which is used in the 
derivation of the gap equation, is that the Schwinger phases of the quasiparticle Green's function and its inverse are 
identical. 



Appendix B: Interaction coefficient functions I n ', n (E) 



The interaction coefficient functions I n ^ n (E) entering the gap equations (f20|) - (|23|) are given by the following 
expression, 



In',n(E) 



2C n ' jTl (y)dy 



y^y 



l-£ 2 



where functions £ n / n (y) are defined in Eq. ([lip and 



arctan 



T,y — 



\_E\_ _ 

huj c y 1 + b 2 y + b 3 y 2 ' 

4tt(1 + hvMv) 



Ky/xy\Jl + b 2 y + b 3 y 2 



+ —(E 2 - 1) + E„ J I - El arctan 



1 - E 2 

y 



E u 



(Bl) 

(B2) 
(B3) 
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In the case of model interaction with a static screening the corresponding coefficient functions I n ', n arc energy 
independent. They can be easily obtained from the more general dynamical expression in Eq. (jBip by substituting 
E = and taking 6j = (i = 1,4), i.e., 

c n ,n'(v)dy (B4) 



^v 7 ^ + 47i-r%) 

Now, in the case of a dynamical screening, the interaction coefficient functions I n ^ n {E) become energy dependent. 
Using the definition in Eq. (|B1|) . we can easily tabulate these functions. The resulting dependence can be fitted by 
Padc approximants of order [2/1]: 

In'AE) = — '-■ (B5) 

Such fits are within a few percent of the numerical results. For the functions I n \n(E) with not very large indices, 
n',n < 10, the fits were obtained from the data in the energy range < E < 12huj c , while for larger indices, 
10 < n! ', n < 40, we were also using data in a wider range of energies, < E < 4QTiuj c . A subset of the corresponding 
fit parameters for I n r, n (E) with < n' , n < 10 in the case of B = 2 T and k — 2 is presented in Table IIIII Note that 
we show only the upper triangular part of the corresponding tables because all coefficients arc symmetric with respect 
to the exchange n' and n. 



Appendix C: Free energy density 

The energy of a given solution is defined by the value of the Baym-Kadanoff effective action calculated on this 
solution. In our analysis, we use the two-loop effective action given in Ref. l3lT Then for the energy of the system, we 
have 



V 



dw , 
2^ 



Tr < — lo 



dG^juj) 



GH + l[rV)GH-i] 



(Cl) 



By making use of the relation 



dG- 1 {u\T,r') 
dw 



5{v-v>), 



we obtain the energy expressed through the translation invariant part of quasiparticle Green's function 

1 



f — tr 

2tt 



-wG(w;0) + - / d?rS- 1 (u);r)G(u);-r) 



-v 0: 



(C2) 



(C3) 



where the overall factor Q is the space volume. Dividing V by this volume and performing all integrations and traces, 
we derive the following expression for the energy density of the system: 



V 



4n£ 2 / 2tt 



9^ 



w - Ms - £A , w - /x s - £A 



(a; - /Xg)(cj + /xgna) + £A Ag ras + - 1) 

(u, + ^„ 5 )2_M|„ s 



V 



i E I 5 - ^ - ^ A °) si § n + 5 ( E hs " A. - ?Ao) sign 

/ — i \ • / v/,/ 2 , r 2 \ , A €n S (A fns +^A )+2a;2 n ( n _i) 



E 

n=2 

E 2w c \]n(n - 1 



(C4) 





1 

2 
3 
4 
5 
6 
7 
8 
9 
10 


1 

2 
3 
4 
5 
6 
7 
8 
9 

10 


1 
2 
3 
4 
5 
6 
7 
8 
9 
10 



1 

2 
3 
4 
5 
6 
7 
<S 
9 
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TABLE III: Coefficients a n n /, 6 n>n < , c n n / and d n n /, evaluated for B = 2 T and k = 2. 



ri = ti' 



n' = 2 


n' = 3 


n'=4 


n' = 5 


n' = 6 


ri = 7 


ri = 8 


ri = 9 


ri 


= 10 


0.2215 


0.1943 


0.1806 


0.1733 


0.1694 


0.1674 


0.1664 


0.166 





1659 


0.2888 


0.2214 


0.1987 


0.1866 


0.1789 


0.1737 


0.1703 


0.1681 





1667 


1.0151 


0.2905 


0.2231 


0.2004 


0.1887 


0.1813 


0.176 


0.1722 





1694 




0.9712 


0.2915 


0.2241 


0.2012 


0.1896 


0.1823 


0.1771 





1731 






0.9381 


0.2918 


0.2246 


0.2016 


0.1899 


0.1826 


o 


1775 








0.9116 


0.2917 


0.2247 


0.2016 


0.1898 





1826 










0.8893 


0.2913 


0.2245 


0.2014 





1896 












0.87 


0.2907 


0.2243 





2011 














0.8532 


0.2901 





2239 
















0.8382 





2894 





















8247 



1.2315 0.2959 
1.0827 



ri = 


ri = 1 


ri = 2 


ri = 3 


ri = 4 


ri = 5 


ri = 6 


n' = 7 


ri = 8 


ri = 9 


ri = 10 


1.186 


0.2911 


0.1683 


0.1186 


0.0915 


0.0745 


0.0631 


0.0548 


0.0486 


0.0437 


0.0396 




0.957 


0.2684 


0.1617 


0.1175 


0.0923 


0.0757 


0.064 


0.0555 


0.049 


0.044 






0.8567 


0.257 


0.1572 


0.1157 


0.0919 


0.0761 


0.0647 


0.0561 


0.0496 








0.794 


0.2491 


0.1538 


0.114 


0.0913 


0.0761 


0.0651 


0.0567 










0.749 


0.2429 


0.1512 


0.1126 


0.0906 


0.076 


0.0652 












0.7143 


0.2379 
0.6862 


0.149 
0.2337 
0.6627 


0.1114 
0.1471 
0.23 
0.6425 


0.09 
0.1103 
0.1455 
0.2268 
0.6249 


0.0757 
0.0894 
0.1094 
0.144 
0.2239 
0.6093 



ri = 


ri = 1 


ri = 2 


ri = 3 


ri =4 


ri = 5 


ri = 6 


ri = 7 


ri = 8 


ri = 9 


ri = 10 


0.0177 


0.0049 


0.0027 


0.0018 


0.0013 


0.001 


0.0008 


0.0006 


0.0005 


0.0004 


0.0004 




0.0138 


0.0047 


0.0027 


0.0019 


0.0014 


0.0011 


0.0008 


0.0007 


0.0006 


0.0005 






0.012 


0.0046 


0.0027 


0.0019 


0.0014 


0.0011 


0.0009 


0.0007 


0.0006 








0.0109 


0.0044 


0.0027 


0.0019 


0.0015 


0.0012 


0.0009 


0.0008 










0.0101 


0.0043 


0.0027 


0.0019 


0.0015 


0.0012 


0.001 












0.0095 


0.0042 


0.0027 


0.0019 


0.0015 


0.0012 














0.0089 


0.0041 


0.0026 


0.0019 


0.0015 
















0.0085 


0.004 


0.0026 


0.0019 


















0.0082 


0.0039 


0.0026 



0.0079 0.0038 
0.0076 



ri = 


ri = 1 


ri = 2 


ri = 3 


ri =4 


ri = 5 


ri = 6 


ri = 7 


ri = 8 


ri = 9 


ri = 10 


0.2047 


0.1273 


0.1041 


0.0902 


0.08 


0.0721 


0.0659 


0.0608 


0.0565 


0.0529 


0.0497 




0.2207 


0.1373 


0.1118 


0.0973 


0.0867 


0.0784 


0.0715 


0.0659 


0.0611 


0.0571 






0.2305 


0.1446 


0.1176 


0.1027 


0.0921 


0.0835 


0.0765 


0.0705 


0.0654 








0.2376 


0.1506 


0.1224 


0.1071 


0.0964 


0.0878 


0.0807 


0.0745 










0.2433 


0.1556 


0.1265 


0.1109 


0.1 


0.0915 


0.0843 












0.248 


0.16 


0.1301 


0.1142 


0.1032 


0.0946 














0.2519 


0.1639 


0.1333 


0.1171 


0.106 
















0.2554 


0.1673 


0.1362 


0.1197 


















0.2583 


0.1703 


0.1388 



0.261 0.1731 
0.2634 
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Note that the sum over Landau levels in the last expression has the following asymptotic behavior at large n: 



1 A 0V V^ X £A^ 



4ir£ 2 



EE 



(C5) 



which is convergent only if the layer- and spin-averaged gap decreases fast enough with n. 

In the framework of the canonical ensemble, the energy of a system is replaced by its free energy through the 
corresponding Lcgendrc transformation. For the system under consideration, we find the following free energy density: 



T = V + nop, 



(C6) 



where 



P = 1 



duj 1 
2?T 



-tr 



G(u/;0) 



Ant 



E 



sign(^ L .J + >'un ;/•;/;.. - 2 ]T sign{fi (ns )0{^ ns - A/f r 



n=2 



(C7) 



Finally, by making use of the explicit expression for the energy density in Eq. (|C4[) . we obtain 



T 



-y\- 



4ir£ 2 



E, 



Cos t Po + sZ - £A ) sign (E, 



L 



) + \ (*%u +P0 + sZ - £A ) sign {E% u 



E 



([i£ na -fJ, - sZ)sign (fj, $ns ) 9 (n\ na - A/| r 



A £ns (£A - A s „ s ) + 2io l c M ins 



M, 



[M^ns — p\ ns ) 



} J 2w cV / n(n - 1) 



(C8) 
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